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Abstract 

The outcome of collisions between small icy bodies, such as Kuiper belt objects, is 
poorly understood and yet a critical component of the evolution of the trans- Neptunian 
region. The expected physical properties of outer solar system materials (high porosity, 
mixed ice-rock composition, and low material strength) pose significant computational 
challenges. We present results from catastrophic small body collisions using a new hy- 
brid hydrocode to iV-body code computational technique. This method allows detailed 
modeling of shock propagation and material modification as well as gravitational reac- 
cumulation. Here, we consider a wide range of material strengths to span the possible 
range of Kuiper belt objects. We find that the shear strength of the target is important 
in determining the collision outcome for 2 to 50-km radius bodies, which are tradition- 
ally thought to be in a pure gravity regime. The catastrophic disruption and dispersal 
criteria, Q* D , can vary by up to a factor of three between strong crystalline and weak ag- 
gregate materials. The material within the largest reaccumulated remnants experiences 
a wide range of shock pressures. The dispersal and reaccumulation process results in 
the material on the surfaces of the largest remnants having experienced a wider range 
of shock pressures compared to material in the interior. Hence, depending on the ini- 
tial structure and composition, the surface materials on large, reaccumulated bodies in 
the outer solar system may exhibit complex spectral and albedo variations. Finally, 
we present revised catastrophic disruption criteria for a range of impact velocities and 
material strengths for outer solar system bodies. 
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Introduction 



Kuiper belt objects (KBOs) are some of the least altered bodies in the solar system. Hence, the 
properties of KBOs and their cometary fragments may faithfully reflect the composition and phys- 
ical structure of the volatile-rich planetesimals that accreted into the larger planets, thus providing 
insights into the conditions in the solar nebula. Astronomical observations of KBOs have revealed 



variations in the surface composition of the largest bodies (Brown 2008). In addition, the colors 



of KBOs may suggest division into subgroups with different dynamical histories (Doressoundiram 



et al.||2001 Jewitt and Luu||2001 Doressoundiram et al.|20 05). Bulk densities and inferences about 
internal structure may also be derived from remote sensing flNoll et aI.|2008||Stansberry et al.|2008j ). 
One of the few modification processes that have acted upon KBOs is mutual collisions, which may 
alter both composition and internal structure (Leinhardt et al. 20081. In this work, we develop 



methods to address the extent to which KBOs have been modified by mutual collisions since the 
last stages of planetary accretion. 

Over the age of the solar system, the surfaces of KBOs have been modified by high energy 
particles, micrometeorite bombardment, impact cratering events and mutual collisions, and possible 



endogenic resurfacing (Stern 2003 1 . While space weathering processes tend to darken and redden 



the surfaces of asteroids (Chapman 2004), the net effect in the Kuiper belt has been difficult to 



quantify, partly because of the unknown surface composition and chemistry. However, the observed 
high albedos and range of colors demonstrate that space weathering does not dominate the surface 
properties of KBOs. Similarly, the role of endogenic processes is poorly constrained. Unlike the 
largest bodies in the asteroid belt, the effects of radioactive decay are minimal or negligible in the 
thermal evolution of Kuiper belt objects because of the longer accretion time scales. As a result, 



few classical KBOs may have differentiated or ever supported cryovolcanic processes (Merk and 



Prialnik 2006). 



The cumulative effects of impact cratering and catastrophic collisions in the Kuiper belt have 
also been difficult to quantify. KBOs currently reside in a low number density region of the solar 
system. So, although the eccentricities and inclinations of the bodies are relatively high in com- 



parison to the major planets, the present day mutual collision probabilities are low (Durda and 



SternpOOO 1 . However, there is abundant evidence for a significant collisional history in the Kuiper 
belt from observations of dust production, the size distribution of bodies, and the low present-day 



total mass of KBOs (Leinhardt et al. 2008). Based on models of the dynamical history of the 



Kuiper belt, most KBOs with diameters < 100 km have suffered large collisions over their lifetime 



(Durda and Stern 2000; Farinella and Davis 1996). It is these large collisions with other KBOs 



that may have significantly changed the internal structure and/or contributed to devolatilization 
of the primordial bodies. In order to assess the effects of collisions on the properties of KBOs, the 
outcome of collisions must be related to observable parameters, e.g., surface composition, colors, 
and densities. The range of possible collision outcomes depends on the initial properties of KBOs 
(mass, density, composition, and internal structure) and the impact conditions (velocity, impact 
angle). Single impact events are generally expected to brighten the surface by excavating through 
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the space-weathered surface into fresh volatile-rich layers. However, if KBOs are highly porous 
and mixtures of ices and refractory materials, small cratering events may result in compaction of 



surface layers and possibly net devolatilization, which could darken the surface (Leinhardt et al 



2008 Housen and Holsapple 2003 ) . The cumulative effects of individual impact cratering events 
will depend on the dynamical history in the Kuiper belt. More energetic impact events lead to 
total disruption of the bodies and gravitational reaccumulation into many smaller fragments. The 
fragments may or may not reflect the general physical properties of the original target. 

The outcome of collision events fall into three general categories: cratering, shattering, and 
dispersal. Cratering events affect only the surface of the target. Shattering events break a coherent 
target into fragments. Dispersal events not only shatter a coherent target, but also cause fragments 
to escape the largest remnant. Each type of collision outcome may produce changes in the physical 
properties of small bodies in the outer solar system ( Leinhardt et al.]|2008 ), but the most energetic 



dispersal events produce potentially observable changes in the internal structure and surfaces of 
the reaccumulated bodies. 

Shattering and dispersal impacts are often discussed with respect to the catastrophic shattering 
criteria and the catastrophic disruption and dispersal criteria, Q* s and Q* D , respectively. The Q* s 
curve describes the energy per unit mass of the target delivered by the projectile such that the largest 
intact fragment is half the mass of the original target. In comparison, the Q* D curve describes the 
specific energy necessary from an impact such that the largest remnant is half the mass - either as 
a single intact fragment or as a reaccumulated rubble pile. In the strength regime (target radii, Rt, 
< 1 km), Q* s and Q* D are equivalent, and Q* s decreases with increasing target size due to the effects 
of strain rate and length scale on the tensile strength ( Grady and Kipp|1980 Housen and Holsapple 
19991. The slope of the Q* s curve changes sign when the lithostatic pressure is comparable 



1990 



to the tensile srength of the target (Rt > 10's to 100 km) (e.g., Melosh and Ryan||1997 ). However 



the slope of the Q* D curve changes sign at a smaller target radii (~ 1 km), when gravitational 
reaccumulation becomes important in determining the mass of the largest post-collision remnant 
defining the transition from the strength to the gravity regime. Here, we focus on the gravity 
regime. The intercept between the two regimes defines the weakest material in the system, which 
is important for the collisional evolution of a population. 

To accurately simulate collisions between Kuiper belt objects, an ideal numerical technique 
must include the capability (i) to use multiple equations of state in compositional mixtures, (ii) 
to model compaction of pore spaces, and (iii) to model a variety of internal structures. In this 
work, we develop a new technique that can more accurately model the possible range of material 
properties of KBOs than techniques used in the past. In this paper, we focus on the effects of 
material strength on Q* D , leaving the investigation of porosity and mixed materials for future work. 
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Numerical Method 



Most of the numerical simulations on catastrophic or near catastrophic collisions have focused 
on the asteroid belt (for example, Benz and Asphaug 1999; Melosh and Ryan]|1997 ) and the special 



case of the formation of Pluto and Charon (Canup 20051. As a result, we will begin our study 



by comparing the results from our hybridized method to the catastrophic collisions of asteroid-like 
rocky bodies before investigating more complex icy physical structures. 

In high velocity collisions, the shock deformation and the gravitational reaccumulation phase 
have very different dynamical times. As a result, it is difficult to model the entire problem with 
one numerical method. In the first phase, the time scale for the propagation of the shock wave is 
determined by the target size divided by the sound speed of the material (few to 10s of seconds for 
a nonporous rocky target with Rt < 50 km). In the second phase, the time scale for gravitational 
reaccumulation is proportional to 1/y/Gp, where G is the gravitational constant and p is the bulk 
density (hours for p = 1 to 3 g cm -3 ). In order to span the different time scales in our investigation 
of catastrophic disruption, we use a hybridized method that combines a hydrocode to model the 
shock deformation and an iV-body gravity code to model the gravitational reaccumulation. 

Hybridized numerical methods have been used by several groups to investigate the outcome 



of high velocity impacts in the asteroid belt over very long time scales. Michel et al. (2001 2002 



2003 2004) and Nesvorny et al. (2006) used a smoothed particle hydrocode (SPH) coupled with 



an iV-body gravity code to model asteroid family- forming events. Durda et al. (2004) used the 



method to investigate binary asteroid formation. All of these groups used effectively identical 



numerical methods which includes using the same SPH (Benz and Asphaug 1995) and iV-body 



codes (Richardson et al. 2000). The SPH code is used to model the impact, the propagation of 



the shock wave, damage accumulation (a proxy for fracturing), and phase changes. Once the shock 
wave had decayed and damage accumulation was complete, the simulation was handed off to the 
iV-body component to model the gravitational reaccumulation phase by directly translating each 
SPH particle into an iV-body particle. 

In this study, we model the shock deformation with an Eulerian shock physics code, CTH 
( McGlaun et al.|1990 ), instead of the Lagrangian SPH code used in previous work. Several features 
of the CTH code are preferred for an investigation of outer solar system objects: an adaptive mesh 
optimizes the resolution of the problem, mixed cell thermodynamics allows for multiple materials 
in each computational cell, and pore compaction models account for the thermodynamics of high 
porosity materials. The gravitational reaccumulation phase of the collision is modeled using the 



iV-body gravity code, pkdgrav (Richardson et al. 2000 Leinhardt et al. 2000 Stadel 2001), the 



same code used in previous hybrid studies. The combination of a shock physics code and an iV-body 
code allow detailed modeling of the impact as well as gravitational reaccumulation. The numerical 
technique is summarized in the example calculation shown in Fig. [T] 
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2.1. CTH 

CTH is a widely-used multi-dimensional and multi-material shock wave physics software pack- 



age developed by Sandia National Laboratories (McGlaun et al. 1990). A two-step solution is 
implemented to move material through an Eulerian mesh. First, a Lagrangian step allows the 
computational cells to distort to follow material motion. During the Lagrangian step, finite volume 
approximations of the conservation equations (mass, energy, and momentum), the material equa- 
tions of state, and the constitutive equations are solved simultaneously. The velocities, energies, 
stresses and strains are updated at the end of this phase. Second, a remesh step maps the distorted 
cells back to the Eulerian mesh and a second order convection scheme is used to flux all quanti- 
ties between cells. In this paper, we use two versions of CTH (7.1 and 8.0); both versions produce 
similar results. Simulations were carried out in either 2D cylindrical symmetry or 3D rectangular 
geometry. 



2.1.1. Adaptive Mesh 

In order to increase the efficiency of the code, CTH also includes an adaptive mesh refinement 
(AMR) feature that changes the computational mesh at user-specified times and locations according 
to user-defined criteria ( Crawford|1999 1 . For example, the primary criteria are designed to increase 



the mesh resolution at the moving locations of (i) free surfaces and material interfaces and (ii) 
steep pressure gradients. The mesh is refined using an adaptive number of blocks in the problem 
domain. Each block is comprised of a fixed number of computational cells. In this work, we 
use 10 cells in each dimension for 2D and 8 cells in each dimension for 3D simulations. In all 
simulations, the smallest block size is determined by requiring at least 20 cells across the radius 



of the projectile in order to sufficiently resolve the peak shock pressures (Pierazzo et al. 1997 1 . 
Refinement subdivides the length scales of a block by factors of 2. Either 2 or 3 levels of refinement 
are used in each simulation depending on the differences in scale between the zeroth level block size 
and the smallest block size. This method allows the resolution to be focused in the most dynamic 
area of the collision and at material interfaces (including void spaces). Fig. [TJ row one, shows the 
evolution of the adaptive mesh for an impact between two basalt spheres at 1.8 km s _1 . In this 
example, the mesh is refined in front of the shock wave, with moving free surfaces, and over the 
entire surface of the projectile. Note that each box drawn in the figure represents one block. 



2.1.2. Equations of State 



In this work, the targets and projectiles are composed of either nonporous basalt or water ice. 
In both cases, we use a SESAME-style tabular equation of state (EOS) flHolian|1984l ). In a SESAME 
table, the internal energy and pressure are tabulated over a density and temperature grid. The 
tabular format has the advantage of explicitly representing complex phase diagrams that cannot 
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be captured with standard analytic models and the disadvantages of being limited by the finite 
grid resolution and the developer-defined range of applicability (e.g., a table may only represent 
the EOS for a single phase). The basalt table is a density-scaled version (p = 2.82 g cm" 3 ) of the 
SESAME table for a-quartz, which includes the phase transformation to stishovite, melting and 
vaporization with dissociation. The table is part of the CTH package and the development of the 



table is described in Kerley (19991. For simulations of impacts into H2O ice, however, pre-existing 
equations of state models and tables were inadequate for our needs (e.g., did not include more than 
one solid phase or did not accurately represent the melting curve). Hence, we developed a new 
SESAME table for H2O that includes vapor, liquid, and three solid phases (ices Ih, VI, and VII) 
that is based on experimentally determined phase boundaries, shock Hugoniots, and thermodynamic 



properties (Wagner and Pruss|2002 Frank et al.|2004 


Stewart and Ahrens|2005 Feistel and Wagner 


2006 


). For a full description of the 5-phase H2O EOS, see Senft and Stewart 


(20081. Because the 



thermodynamic properties of H2O are difficult to represent with standard analytic formulations 
(such as Tillotson ( |Tillotson||l96"2"l ) and ANEOS ( |Thompson and Lauson|l972l )), we decided to use 
published equations of state for the individual vapor, liquid and solid phases to build an equation 
of state table that is applicable to impact events over the entire solar system. Note that our H2O 
equation of state will produce much more reliable results compared to previous work using the 
Tillotson or ANEOS equation of state because it includes the crucial solid phase transformations 



to ices VI and VII that control the criteria for melting ice with a shock wave (Stewart and Ahrens 
20051). 



2.1.3. Strength Parameters 



The failure of brittle materials (e.g., rocks and ice) is well described by a pressure-dependent 



yield surface (e.g., Jaeger and Cook 1979), where the shear strength has a value of Yq at zero 
pressure and increases to a maximum value of Ym with increasing pressure with a slope defined by 
the coefficient of internal friction (pi. Fragmented materials have lower zero-pressure shear strength 
(cohesion, Y c ) and friction coefficient (4>d) compared to intact materials and the same limiting Ym- 

In all except the hydrodynamic simulations, we use a pressure-dependent shear strength model, 



either the geological yield surface (GEO model, Crawford et al. 2007) or the ROCK model (Senft 



and Stewart 2007). The yield surface describes the amount of shear stress that may be supported 



by the material as a function of ambient pressure. When a measure of the shear stress (the root 
of the second invariant of the deviatoric stress tensor, y^) is greater than the shear strength, the 
material is failing in shear and undergoes irreversible plastic deformation. At lower stresses, the 
material responds elastically. Note that under planar shock wave conditions with uniaxial strain, 
the maximum shear stress, (o"i — 0"3)/2 (where a is a principal stress component), is equal to 

In the ROCK model, shear strength is linearly degraded from an initial intact curve to a 
fragmented curve (friction law) using a dimensionless scalar measure of damage. Damage represents 
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the fraction of integrated effective plastic strain to failure, which is a proxy for the degree of 
fragmentation; hence, intact rock has zero damage and completely fragmented rock has a damage 
of one. The effects of damage are neglected in the GEO model. In both strength models, the shear 
strength is degraded to zero as the temperature approaches the melting temperature (T m ). 

When negative stresses exceed the fracture (tensile) strength, the material fails in tension 
and cracks. In CTH, tensile failure is modeled by inserting void space into a computational cell 
until the volume-averaged stress equals the tensile strength, Yt- In the ROCK model, damage 
also accumulates with tensile fracturing via a crack propagation model. Note that both shear and 
tensile failure lead to fracturing in brittle materials. 

Because of the large uncertainty in the physical properties of planetesimals, we consider a wide 
range of strength parameters for materials with no significant porosity (note that natural intact 
rocks generally contain several percent porosity) (Table [T]). A large range of physical properties 



is expected based on the inferred range of bulk densities of outer solar system bodies (Noll et al 



2008). The parameters chosen are meant to represent general classes of materials, rather than a 
specific rock type. For example, a monolithic body may have similar bulk density and average 
EOS as a densely packed breccia, but the two would have very different shear and tensile strengths. 
A monolithic rock target is represented by the strong strength parameters. A strong model with 
high tension test case is also considered (strongHT). A breccia with a weak matrix is represented 
by the weak strength parameters. An intermediate case (medium) is considered to understand 
the systematic dependence on strength parameters and is not meant to correlate with a particular 
material group. Finally, as an end member, we consider a hydrodynamic (no shear strength) 
material (e.g., a liquid body) with tensile strength. 

The ROCK model parameters are derived directly from laboratory quasi-static and dynamic 



strength measurements (see Senft and Stewart (2007) for basalt and Senft and Stewart (2008) for 
ice). The ROCK model successfully reproduces crater sizes and damage/fracture patterns observed 
in laboratory impact cratering experiments in strong rock ( Senft and Stewart||2007 ). The dynamic 
shear strength of a material is known as the Hugoniot Elastic Limit (HEL), which constrains the 
value of Ym- The HEL of strong rocks is typically a few GPa and the coefficient of friction is 



about one (Melosh 1989). The strong model parameters are based on intact basalt data and are 



similar to values for other strong rocks such as granite (see discussion and references in Senft 



and Stewart 2007 Collins et al. 2004). For comparison, the basalt targets in Benz and Asphaug 



(1999) were modeled with a Von Mises shear strength (no pressure dependence) of Ym = 3.5 GPa. 



The weak model parameters fall in the range observed for terrestrial breccias and cements, with 



HEL values between negligible (unmeasurable) to < 1 GPa (Willmott and Proud 2007 Tsembelis 



et al. 2000 2002, 2004). Note that the HEL of aggregates with a weak matrix is controlled by 



the strength of the matrix (Tsembelis and Proud 2006). The weak model parameters are slightly 



lower than quasi-static data on cold ice and comparable to warm ice dSenft and Stewart||2008[ ); the 



temperature-dependent HEL of H2O ice ranges between 0.05 and 0.6 GPa (Stewart and Ahrens 



2005). The ice targets in Benz and Asphaug (1999) were modeled with a constant shear strength 
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0.1 GPa, purham et al.(l983| ) (see 



of Ym = 1 GPa, much larger than found experimentally (Ym 
also §0. 

Tensile strength has a strong scale-dependence due to absolute length (larger bodies have 
longer fractures; longer fractures are weaker) and strain rate (impacts by larger bodies have lower 
strain rates; brittle materials fail in tension at lower stresses at lower strain rates). As a result, 
the tensile strength decreases with increasing size. We extrapolate laboratory-scale values for 



-(0.25-0.33) 



tensile strength to 1-10 k m-sized bodies following the method of Housen and Holsapple (19901 
where Yr ~ R 



(Housen and Holsapple 1999 Ryan and Melosh 1998| . Notethat the 



uncertainty in the slope of the size dependence leads to an order of magnitude range of values for 
tensile strength at km scales (and a larger range for larger sizes) (see Housen and Holsapple| [l999 ) . 
We use values for the tensile strength that falls within the predicted range for weak and strong 
materials. Typical laboratory-scale dynamic tensile strengths for strong rock are order 10 8 Pa (e.g., 



Cohn and Ahrens 1981 1 and decrease to order 10 6 to 10 7 Pa for 1-10 km sized bodies. The strongHT 



model, with a laboratory-scale tensile strength, tests the sensitivity of the disruption calculations 
to variations in tensile strength. The tensile strength of ice is about an order of magnitude smaller 



than strong rock (Lange and Ahrens 1983). 



In the ROCK model simulations, damage initially accumulates primarily via shear failure 
during the propagation of the primary shock wave. At later times, upon reflection of the shock 
wave from free surfaces, damage also accumulates in tension. The total damage variable is the 
sum of the shear and tensile damage ( Senft and Stewart]|2007 l. In all cases except for the smallest 
targets with the highest tensile strength (strongHT model), the largest intact fragment is smaller 
than a computational block (j ]3.1| see also Benz and Asphaug|1999 Ryan and Melosh|1998 ). Hence, 
each computational block represents a size distribution of fragments on similar trajectories. 

Although the ROCK model is a more accurate description of the strength of rocks, the GEO 
model, with far fewer parameters (6 vs. 13), is easier to intepret. In the GEO model, all simulations 
use a friction coefficient of 1, Poisson's ratio of 0.3, and melting temperature of 1474 K. Simulations 
using the ROCK model are labeled as such Tables [T] &; |2j otherwise the GEO model is used. As 
described below, both material models produce quantitatively similar results for gravity-regime 
catastrophic disruption calculations. 



2.1.4- Self- gravity 

CTH also has the capability to include self-gravitational forces in the calculation. The parallel 
tree method of Barnes and Hut (1986) is included in the version 8 distribution. In previous work 
where the intact fragment size distribution is derived, self gravity is important during the shock 
deformation phase for bodies with radii greater than several 10's km, when central lithostatic 
pressure becomes comparable to the tensile strength. For calculations that do not include self- 
gravity, the model tensile strength must be appropriately increased to mimic the effect of lithostatic 
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pressure (Ryan and Melosh 1998 Benz and Asphaug 19991. In the simulations presented here, 
self-gravity has a negligible effect on the velocity evolution during the short CTH portion of the 
calculation. Hence, it is not utilized for most of the calculations because of the computational 
expense. The insensitivity of the results to self-gravity in the first few seconds is verified by a few 
test cases including self-gravity simulations (see 



3.1.2 and sims. 19b and 19d)). 



2.2. pkdgrav 

Once the initial shock wave and the subsequent rarefaction waves have passed through the 
target, it is no longer necessary or practical to continue the simulation with CTH. Almost all of 
the physics can be captured with an iV-body gravity code. In the example shown in Fig. [TJ row 
2 shows the pressure gradient in grey scale versus time. By 60 seconds most of the shock-induced 
deformation is complete and gravity is the dominant force. 

In these simulations, we have employed the parallelized, Lagrangian, hierarchical tree code 



pkdgrav (Richardson et al. 2000 Leinhardt et al. 2000) to complete the gravitational reaccumu- 



lation phase of the collision. Pkdgrav is a second order symplectic integrator using the leap-frog 
integration scheme. The gravitational evolution of particles is modeled under the constraints of 
self-gravity and physical collisions. The particles have a mass and radius, but no explicit equation 
of state. Material properties are captured by the coefficient of restitution. 

The last output of CTH is run through a translator to convert the Eulerian grid data into 
Lagrangian particles, creating initial conditions for pkdgrav (Fig. [T] rows 3 and 4). If the CTH 
calculation is in 2D, the last time step is rotated about the axis of symmetry creating a 3D object. 
Each adaptive mesh block is then an annulus of material. The annuli are divided into spherical 
particles with a radius determined by the smallest of the 2D block dimensions. If the CTH calculation 
is 3D, the adaptive mesh blocks are translated into spherical particles with a radius determined by 
the smallest dimension of the block. To eliminate simultaneous collisions, which pkdgrav cannot 
tolerate, a small amount of random noise is added to the velocity of each particle (1% the escape 
speed) when handing off from CTH. Because each CTH block represents a size distribution of particles, 
no provision is made to bind certain pkdgrav particles together as an intact fragment, although this 
would be necessary to model the catastrophic disruption of bodies of smaller bodies that transition 
from the gravity to the strength-dominated regime. 

If there is only a small fragment of material within a block, a straight translation into a La- 
grangian particle would create a very low density particle with an artificially large radius. Therefore, 
to avoid unrealistic particles, the minimum density of each particle is limited to 0.75 of the original 
bulk density of the target. If the mass of a block is less than 4 x 10 _7 M to t, no particle is created. 
The cutoff value is an empirical threshold that eliminates numerous tiny ejecta particles that do 
not accumulate into the largest remnants. In addition, little ejected material leaves the CTH mesh 
before handoff. Therefore, the total amount of mass discarded in the translation is negligible and 
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usually involves material that is escaping from the largest remnant of the collision event. 



To facilitate translation into spherical particles, the geometry of the adaptive mesh blocks 
are kept close to square or cubical. The total mass of the square/cubic block is conserved upon 
translation, thus the mass density of the pkdgrav particle is about a factor of two higher than the 
mass density in the original CTH block. This means that some assumption must be made about 
the macroporosity of the post-collision remnants in order to determine a radius from the mass of 
the remnant. Note that in the work presented here, negligible amounts of vaporized material is 
produced. Thus, the effects of volume changes are ignored in the translation. 

Each pkdgrav particle is modeled as an indestructible sphere. Collisions between particles 
follow one of three scenarios: perfect merging, conditional merging, or inelastic bouncing. Tangen- 
tial and normal coefficients of restitution are used to govern the amount of kinetic energy lost as 



a result of each particle collision (Richardson 1994 Leinhardt et al. 2000). The impact velocity 
is given by v = v„ + V(, where v n is the component of the velocity vector that is normal to the 
surface of the target particle and is the tangential component. The post-impact velocity is given 
by v' = — e n v„, + e^v^, where e n is the normal coefficient of restitution and et is the tangential co- 
efficient of restitution. Both e n and et have values between and 1. In this paper, perfect merging 
(e n = 0) is used in most cases to reduce calculation time (denoted by 'merge' in the simulation 
type in Table [2]); however, information is lost about the shape and structure of the reaccumulated 
body. In order to examine the structure of a remnant, inelastic collisions (e n is a nonzero constant 
based on the type of material) are used in some 3D cases (see § 3.3). Conditional merging (e„ = 0, 
if the collision speed is less than the escape speed, otherwise e n > 0) was used to verify the results 
of simulations using perfect merging. Recent field observations and friction experiments on rocky 
materials constrain the value of e n ( Chau et al.|2002 ). In general, e n varies between 0.1 and 0.5 for 
rocks, soils, and ice ( |Chau et al. 2002 Higa et al. 1996 1998). In this work, simulations denoted 
by 'bounce' in Table [2] use e n = 0.5 and simulations denoted by 'loweps' use e n — 0.1. In all cases, 
et is fixed equal to 1, which implies negligible surface friction as found in experiments on rocky 



materials (Chau et al. 2002). 



The pkdgrav phase of the calculation continues until the largest post-collision remnant reaches 
dynamical equilibrium. Dynamical equilibrium is defined by the criteria that the total mass ac- 
creting and orbiting the largest post- collision remnant must be < 10% of the mass of the largest 
remnant (Mi r ); in addition, a minimum of 50 dynamical times is imposed upon all calculations. 

In this work, the number of pkdgrav particles ranges from several thousand to several tens of 
thousands (Ni n it in Table [2|. Recall that each pkdgrav particle represents a CTH computational 
block with 8 3 cells in 3D and 10 2 cells in 2D. The AMR feature allows the shock deformation to be 
very well resolved while minimizing the number of particles in the pkdgrav phase of the calculation. 
Hence, the CTH-pkdgrav technique highly resolves the energy and momentum coupling to the target 
at relatively low computational expense. 
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2.3. Tracers 

To track the motion and physical properties of material flowing through the mesh, massless 
Lagrangian tracer particles are included in the CTH calculation to provide material histories, includ- 
ing trajectories, pressures, temperatures, etc. In the 2D CTH calculations, 900 tracer particles are 
placed in the target on rays initiating from the impact site (Fig. (2^). In the 3D CTH calculations, 
tracer particles are placed in a uniform 3D grid in both the projectile and target. At handoff from 
the hydrocode to the gravity code, the tracers associated with each pkdgrav particle are recorded 
for analyses of the material deformation history. In Fig. [2^3, particles that contained no tracer par- 
ticles are assigned a peak pressure history value by averaging the peak pressure from the nearest 2 
neighbors. 

It is difficult to ensure that the divergent flow in the ejecta cone is filled with tracer particles. 
Because it is highly resolved and comprises a significant surface area, the ejecta cone requires the 
most interpolation. However, most of the material in the ejecta cone escapes after the impact. As 
a result, interpolation of a large amount of the ejecta cone will not effect any conclusions that are 
drawn about the properties of the largest remnant ( ^3.3| ). The original surface of the target near 
the impact site also needs to be interpolated because it has been highly resolved by the adaptive 
mesh refinement, resulting in pkdgrav particles between tracer particles. 



2.4. Handoff Time 

In the hybridized method presented here, the transition from hydrocode to iV-body gravity 
code occurs after the shock and rarefaction wave has propagated through the system and dissipated. 
After this point, modification of material due to shock waves is negligible. Figure [3] shows the mass 
of the largest post-collision remnant versus time of handoff from CTH to pkdgrav for an impact 
between a 10 km radius target and a 0.83 km radius projectile at 1.9 km s" 1 (sim 8c in Table [2|. 
The mass of the largest remnant levels off between 10 and 15 seconds or about four to six times the 
sound crossing time. A least squares fit of the data after 10 seconds finds a slope consistent with 
zero (—0.002 ± 0.003). There is scatter in the mass of the largest remnant due to varying particle 
number and particle size from simulation to simulation. The resolution of the iV-body component 
of the simulation varies depending on the handoff time since CTH utilized adaptive mesh refinement. 
We estimate a conservative error in the mass of the largest remnant to be ±0.1 M; r /My, where 
Mt is the target mass, based on the peak to peak scatter of the test results. In the rest of the 
simulations presented in this paper, the handoff time was > 4 times the sound crossing time of the 
target. 
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Results 



We conducted sets of simulations with (i) fixed projectile-target size ratios and varying impact 
velocity or (ii) fixed impact velocity and varying projectile size. In this work, we are primarily con- 
cerned with the mass and physical properties of the largest remnant. The critical disruption energy, 
Q* D , is determined by a non-linear least squares fit to the resolved largest remnants (M^/Mt > 0.1). 
Table [2] presents a summary of the simulation parameters and results. The simulations denoted 
with an asterisk are close to or below the resolution limit of the hybrid calculations (0.1 M^/Mt) 
and have a low number of pkdgrav particles (Ni r < 200). (If it is necessary to use an under-resolved 
simulation to determine Q* D the under-resolved simulation is given less weight than the resolved 
simulations). Tables [3] & [4] present the critical impact velocity (for a fixed projectile to target mass 
ratio) and critical projectile radius (for a fixed impact velocity) and the corresponding Q* D values 
for each set of simulations (Fig. Q. We vary the strength, impact angle, pkdgrav collision model, 
and composition to determine their effects on the energy coupling in catastophic disruption events 
and the physical properties of the largest reaccumulated remnants. 



Next, we describe the validation of the hybrid method by comparison to previous work (§ 3.1 ). 
We conduct simulations with parameters comparable to results on the catastrophic disruption of 
strong bodies presented in Melosh and Ryan (1997) (sims. 2, 10, & 16) and Benz and Asphaug 



(1999) (sims. 7 & 12). In section § 3.2 we describe the results of varying the shear and tensile 



strength of the bodies while holding the target to projectile size ratio fixed (Table [l}. In 3D 
simulations, we consider the effects of impact angle and composition on the physical properties of 



the largest remnant (§ 3.3). Finally, we examine the role of the pkdgrav collision model on the 



gravitational reaccumulation process (§ 3.4). 



3.1. Method Validation 



The catastrophic disruption criteria derived from our hybrid calculations are presented in 
Figure U and Tables |] & g). In Fi gure El the results are compared to the 90°, strong target 
simulations by Benz and Asphaug (1999) (solid lines) and Melosh and Ryan ( 1997[ ) (dashed lines). 
Although the details of the shear and tensile strength models are different, the strong target results 
(green and red symbols) are in good agreement with previous work. The weak and hydrodynamic 
targets (black and blue symbols) have disruption criteria that are significantly smaller than the 
strong target cases (see 



3.2). 



In this work, we held the size ratio between the target and the projectile fixed for most 
simulations rather than holding the impact velocity constant in order to ensure Mp <C My. With 
the fixed target and projectile sizes, the optimum resolution for the calculations can be easily 



controlled. In previous work with a fixed impact velocity (e.g., Benz and Asphaug 1999), the size 



and resolution of the projectile varied significantly (and in some cases Mp > Mt)- In simulations 
2, 10, and 16 (red points in Fig. HI), the size ratio between the projectile and targets are the same 
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as in Melosh and Ryan (1997 1 . The derived critical impact velocities (that correspond to the fitted 
Q* D in Table 3) are 1.4, 4.8, and 1.9 km s _1 , respectively, compared to 1.3, 3.7, and 1.8 km s~ ] 



m 



Melosh and Ryan ( 1997 ) . Although our results are systematically higher, perhaps indicating that 



our shear strength parameters are slightly stronger than the unreported values used by Melosh and 



Ryan (19971, the critical impact velocities are generally within our error bars. Simulations 7 and 
12 have fixed impact velocities of 3 km s _1 , as in Benz and Asphaug (19991, and our results agree 



within error (open green squares and solid line in Fig. |4J) . Note that there is comparable scatter in 
the individual Q* D data points that are used to derive Q* D curves between our and previous work. 

Another measure of the accuracy of our hybrid simulations is the mass of the largest remnant 



as a function of impact energy. Benz and Asphaug (19991 found that the normalized mass of the 



largest remnant, Mi r /Mx, follows a tight linear trend with the normalized impact energy, Q/Q* D , 
in the hypervelocity regime. For example, they fit a slope of —0.5 for 3 km s _1 impacts into basalt 
targets. We find a similar slope of —0.48 ± 0.02 (Fig. [5]) for impact velocities over the larger range 
of 0.7 to 5.6 km s _1 into basalt targets. The linearity of the results justify the use of a linear fit 
to the remnants from each simulation within a group to derive the value of Q* D . Note that the 
slope would be -1 for pure energy scaling (Stewart & Leinhardt, in prep.). The slope and linearity 
of the largest remnant mass is a satisfying validation of the hybrid technique. The tight linear 
relationship is independent of the impact velocity and the strength of the targets, which varies 



from hydrodynamic to strong. Benz and Asphaug (19991 noted a slight dependence of the linear 



slope on impact velocity, which may include differences in the projectile to target size ratio. Benz 



and Asphaug (19991 also observed slightly steeper linear slopes for the normalized mass of the 



largest remnants in ice targets. 

Hence, we confirm that the CTH-pkgrav hybrid technique is in excellent agreement with pre- 
vious studies of catastrophic disruption in the gravity regime. 



3.1.1. Comparing Different Impact Velocities 



It is well established that changing either the impact velocity or the mass ratio between the 
projectile and target will result in differences in the energy and momentum coupling to the target 
(e.g., Housen and Holsapple 1990 Melosh and Ryan 1997 1 . Therefore, detailed comparison of 



catastrophic disruption results between different studies and between different target sizes requires 
an explicit correction for the different modeled impact velocities and/or mass ratio. In order to 
be able to extend our results to different (hypervelocity) impact velocities (where Mp <C Mt), 
we apply a impact velocity correction using the fragmentation theory developed by |Housen and| 



Holsapple (1990). 



First, we use the slope of the Q* D curve to constrain a material parameter, fj,, in the coupling 
theory of Housen and Holsapple (1990). In the gravity regime, Q* D would be proportional to Bjj> if 



material properties were inconsequential and pure energy scaling applied. In reality, the coupling 
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of the projectile's kinetic energy to the target depends on the material and the impact velocity. In 



the analytical model of Housen and Holsapple (1990, Eq. 36) for the critical shattering criteria, Q* s , 
material properties are captured by a coupling parameter that scales with rather than scaling 
with the impact energy. The exponent, fi, is a simple descriptor of material properties with a value 
of 0.55 to 0.60 for rocks and water and ~ 0.4 for sand ( Housen and Holsapple|1990 Holsapple|1993 



Schmidt and Housen 



(Housen and Holsapple 1990 Eq. 68). 



19871). In the gravity regime, Q* s ~ R^ ~ R 



,1.2-1., 



for sand to rocky targets 



Next, we assume that the additional energy required for dispersal of the shattered fragments 
produces a negligble change in slope between the Q* D and Q* s curves. This assumption is supported 



by Housen and Holsapple (1990 Eq. 73) 's estimate of Q* D based on the fraction of shattered 
fragments that reach escape velocity. Using laboratory data as a constraint on the fragment velocity 
disruption, the exponent of the Q* D curve is negligibly smaller than the Q* s curve in the gravity 
regime. Numerical simulations of catastrophic disruption support the departure from energy scaling 



and similarity in slope to Q* s . The 90°, basalt, gravity regime data from Benz and Asphaug (1999) 
(Vi =3 km s _1 ) and Melosh and Ryan (1997) (Vi =1.3-5.3 km s _1 ) have fitted exponents of 1.26 



and 1.44, respectively, implying that the value of fi falls in the range of 0.42 to 0.48. 

Laboratory and numerical experiments demonstrate that increasing the impact velocity in- 



creases the critical disruption energy for a fixed target radius (Housen and Holsapple 1990 Benz 



and Asphaug][l999l ). With increasing hypervelocity impact velocities, more of the impact energy is 



expended in material deformation and the disruption process becomes less efficient. Housen and 
Holsapple ( 1990 Eq. 68) find that Q* s scales by yr 3 ^+ 2 in the gravity regime, and we assume that 



the same velocity correction applies to Q* D . The derived correction factor is in reasonable agree- 



ment with results from Benz and Asphaug (1999), Q* D increases by a factor of ~1.7 between lines of 
constant Vi = 3 km s _1 and Vi = 5 km s _1 for basalt targets. In addition, in our own simulations, 
the velocity corrected fixed mass ratio simulations 6 & 11 are consistent with the constant impact 
speed varying mass ratio simulations of the same strength model 7 & 12. The original critical 
impact velocities (pre-correction) are noted in Table [3j 

We illustrate the velocity correction in Figure|4|B, where the Q* D data are adjusted to a constant 
impact velocity of 3 km s -1 . The correction uses an empirical value of \i ~ 0.45 based on a mean 
gravity regime slope (3/j, ~ 1.35) from our and previous worker's results. Our hybrid simulations 
of strong and weak basalt targets are fit with slopes of 1.3 and 1.2 respectively (dotted lines in 
Fig. |4j3). Interestingly, the two hydrodynamic cases suggest a slope of 1.9, approaching the ideal 
Q* D ~ Rj, for pure energy scaling in the gravity regime. 

Note that the dependence of Q* D on shear strength is confirmed to be a primary result and not 
an artifact of the differences in critical impact velocities. The physical explanation for the role of 
shear strength is dicussed in § |3.2 
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3.1.2. Other Checks 

We confirm that each pkdgrav particle represents a size distribution of smaller fragments. The 
targets and projectiles are completely damaged {D = 1) for all simulations except the strongHT 
cases with Rt = 2 km. In the strongHT case, the tensile strength is the highest expected for 
km-scale bodies. These results are consistent previous work, where the mass fraction of the largest 
intact fragment in the reaccumulated remnant drops precipitously above target radii of a few 
100 m (Benz and Asphaug 1999 Melosh and Ryan|1997 ). In the 2-km strongHT cases (sims. 2, 6R, 



7R), the largest fragment is composed of several pkdgrav particles that all merge into the largest 
gravitationally reaccumulated remnant. Some of the the scatter in the Rt = 2 km results reflects 
the increased influence of shear and tensile strengths compared to the larger bodies. However, the 
amplitude of the strength effects at this size are within error within each strength model group. 

In simulations 1-16, we use 2D geometry in CTH to decrease the computation time. We verify 
that calculations using 2D geometry are robust by comparing identical impact conditions in 2D 
and 3D geometries. The mass of the largest remnant in simulation 17a {Mi r /Mx = 0.57) is nearly 
identical to the mass of the largest remnant in simulation 14c (M^/Mj 1 = 0.60). Therefore, we 
have determined that simulations using 2D and 3D CTH components are in very good agreement 
with one another. 

Finally, we confirmed that self-gravitational forces are negligible during the CTH component 
of the calculation. Simulations 19b and d include self-gravity in the CTH component. The mass 
of the largest remnant and the second largest remnant is virtually unchanged when compared to 
simulations without self gravity in the CTH component (sim 19a and c). In addition, the resolution 
in simulations 19c and d are higher than in simulations 19a and b. Once again, the masses of the 
largest post-collision remnants are indistiguishable. 



3.2. Strength Effects on Q* D 



Small bodies in the solar system possess a wide range of shear strengths, from monolithic 



bodies inferred from fast spin rates, (e.g., Pravec and Harris 20001 to nearly strengthless bodies 
(e.g., ~ 1 — 10 kPa for comet 9P/Tempel 1, Richardson et al. 2007 ) . Previous simulations have 
shown that the catastrophic shattering criteria (Q* s ) is typically two orders of magnitude less than 



the criteria for dispersal (Q* D ) in the gravity regime (Melosh and Ryan 1997 1 . As a result, the 



role of strength has largely been ignored in the gravity region of the Q* D curve. Although previous 
workers have considered a wide range of materials, including ice, basalt and mortar, the literature 
has focused on the differences in tensile strength and not considered a wide range of shear strength 



(Ryan and Melosh 1998 Benz and Asphaug 1999 1 . 



In this work, we consider a wide range of shear and tensile strengths. We find that the 
catastrophic disruption criteria is higher in targets with higher shear strength. As expected, the 
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effect of tensile strength is negligible (Fig. [4|. An increase in the tensile strength by two orders 
of magnitude (red symbols: = 10 5 Pa, green symbols: Yj* = 10 7 Pa) has a small effect on the 
collision outcome at a target radius of 2 km (within error) and essentially no effect on the collision 
outcome at a target radius of 10 km. In contrast, the shear strength model has a large effect on the 
resulting Q* D , with variations up to a factor of 8 at Rt = 10 km between weak and strong targets. 

The shear strength affects the catastrophic disruption criteria in the gravity regime by modi- 
fying the pressure and velocity distributions in the target. Shear strength increases the peak shock 
pressure and shock wave decay rate compared to a strengthless material. In Fig. [G]^, the peak 
shock pressures along a ray of tracer particles adjacent to the centerline (refer to Fig. [2j are shown 
for identical impact conditions but varying shear strength (sims. 8a, 10c, 11c, and 13a). The peak 
pressure is higher with greater shear strength because the shock Hugoniot is shifted upward by 2/3 
of the yield strength compared to a hydrodynamic Hugoniot ( Melosh||1989 ). However, as the shock 



wave propagates into a target with higher shear strength, a portion of the energy is partitioned 
into plastic deformation and the shock wave decays more rapidly compared to a weaker material. 
In other words, the stronger the material, the more energy is partitioned into overcoming the shear 
strength. 



The cumulative distribution of peak pressures reflects the different shock decay profiles (Fig.[6p). 
Because of the shallower shock decay rate in the weak target, a larger mass fraction reaches higher 
shock pressures compared to the strong targets. The fraction of highly shocked material decreases 
as the strength of the material increases. The initial velocity evolution after the impact is con- 
trolled by the shock pressure field. Because the weak target has a shallower shock decay profile, 
the pressures are higher at the time the shock wave encounters the free surface. The material 
velocity follows pressure gradients, so the velocity distribution in the weak target is faster com- 
pared to the strong targets (Fig.[6p). Thus, more material reaches escape velocities and the largest 
reaccumulated remnant is smaller for weak targets. 

Note that the pressure profiles using the GEO model are very similar to the ROCK model 
(e.g., sims. 10c & 11c in Fig. [6j^), which leads to similar Q* D values for each of the strength groups 
(see Table [3] weak - 1, 4R, strong - 2, 5R, strongHT - 3, 6R). As expected, the pressure decay 
profile in the medium strength target falls in between the weak and strong target cases. 



Using a physical model derived from impact cratering, Melosh and Ryan (19971 used the 
pressure decay profile to predict the catastrophic disruption criteria, assuming an average power 
law decay exponent of -2. The material velocity at the antipode from the impact is derived from 
the peak pressure using the conservation equations and assuming velocity doubling at the free 



surface. Melosh and Ryan (1997 1 's catastrophic disruption criteria was assumed to correspond to 
an antipodal velocity of half the escape velocity of the target. From Fig. [6]^, it is clear that shear 
strength has a significant effect on the pressure decay profile and the resulting material velocities. 
The antipodal velocities range from about 20 m s~ x for the strong target cases to 40 m s~ x for 
the weak target. Both cases are greater than the escape velocity, even though the strong target 
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case is subcatastrophic. Although the conceptual model from Melosh and Ryan (1997) is a very 
useful approximation for a disruption criteria, the true velocity distribution is more complicated 
and antipodal velocities alone are insufficient to predict disruption. 

We have found that the shear strength affects the catastrophic disruption criteria by modifying 
the pressure and velocity distributions in the target. In this work, we considered general groups 
of material properties. We plan to conduct a more focused investigation on the specific properties 
of pure ice and ice-rock mixtures in a later paper using ROCK model parameters for ice that are 



currently under development (Senft and Stewart 2008) 



3.3. Properties of the Largest Remnant 

Because the gravitational reaccumulation stage is calculated in our hybrid technique, we are 
able to examine the physical properties of the reaccumulated remnants. Although we are interested 
in all of the collision remnants, it is the largest remnant that has the highest numerical resolution. 
The physical properties of large collision remnants may also be observed astronomically. For exam- 



ple, some of the largest Kuiper belt objects are inferred to be rubble piles (Trilling and Bernstein 



2006 Lacerda and Jewitt 2007). 



First, we consider the range of shock deformation experienced by the material that reaccu- 
mulates into the largest remnant. In Fig. [6p, the peak shock pressure distribution is shown for 
the largest remnant in the four simulations with identical impact conditions but different strength 
models. When compared to the shock pressure distribution of the total mass (Fig. (6^3), we see that 
the largest remnant materials have experienced a large range of shock pressures in the collision 
event. In general, as expected, the shock pressures are lower in the largest remnant, but some 
fraction of the remnant is derived from the most highly shocked material. Although peak shock 
pressures are not particularly high compared to other hypervelocity impacts in the solar system, 
if the material responds in an irreversible way to low shock pressures (e.g., by structural changes, 
phase changes, or chemical changes), some of the materials in the largest remnant will reflect these 
changes. 

Next, we examine the effect of impact angle on the largest remnant. We have conducted most 
of our simulations at 90° for computational efficiency, but in general, an impact between two small 
bodies in the outer solar system will not be head-on. In order to illustrate the effect of impact 
angle, we compare two impact events on a 50-km radius target that have the same mass largest 
post-collision remnant {Mi r /M.T ~ 0.7). Figures [T]^ and B show the cumulative peak pressures 
attained in a 90° (sim. 14d) and a 45° (sim. 18a) impact. The solid line is peak pressure attained 
by the total mass, the dashed line is the largest remnant, the dotted line is the second largest 
remnant. Both the largest and second largest remnants show a similar range in peak pressure 
distribution as the total mass. Note, however, that the second largest remnant is generally under- 
resolved in our calculations. The finding that material in the largest remnant experiences the full 



- 18 - 



or nearly full range of peak shock pressures is a general result from all of our simulations. The 
cumulative pressure distribution is wider, meaning a larger range of shock pressures, for the 45° 
case compared to the 90° case. The 90° impact has a larger mass fraction of material in the high 
pressure tail, while the 45° impact has a larger mass fraction in the low pressure tail. Because an 
oblique impact has less efficient energy coupling to the target, the impact velocity must be faster to 
reach the same outcome (same M^/Mt)- As a result, some of the material reaches a higher peak 
pressure. 

We also examine the provenance of material that reaccumulates into the largest remnant in 
catastrophic-level collisions, where M^/Mt ~ 0.5. Previous studies have noted that most material 
originates from the antipoldal region to the impact, but not examined the source and reaccumulation 
in more detail. Fig. [8] shows the original (pre-collision) location of particles that reaccumulate into 
the largest remnant for 3D simulations of 50-km radius targets (sims. 17-19). Particles are color 
coded for peak shock pressure. The left column shows simulations where the perfect merging 
criteria is used in the pkdgrav phase. The right column shows the same simulations with inelastic 
collisions (bouncing) with a normal coefficient of restitution of 0.5. The striking differences between 



the merging and bouncing simulations are discussed in § 3.4 The true provenance of material in 



the largest remnant should fall between these two limiting cases. 

In all simulations, the largest remnant is primarily composed of material from the antipodal 
hemisphere suggesting any compositional diversity from the original target could be preserved in 
the largest remnants in catastrophic-level collisions. Note that the pressure range in the ice target 
is lower compared to the basalt targets. The lower shock pressure range results in a generally lower 
material velocity distribution. In the gravity regime, ice targets disrupt at lower specific impact 
energies, primarily because of the lower density of the ice and the resulting smaller escape velocity. 

Although the interior structure of small bodies in the outer solar system cannot be observed 
directly, their surface properties can be studied and internal structures can be inferred. The results 
from our simulations can help to determine whether features observed on the surface can be ex- 
tended to the interior. Fig. [9] shows an example of a largest remnant from a 3D, 45° impact into a 
50-km radius basalt target (sim. 17b). The exterior of the remnant is heterogeneous with particles 
that have a broad range of peak pressures (colors ranging from yellow to purple). The interior 
of the largest remnants is slightly less heterogeneous - composed mostly of moderately shocked 
material (mostly green). This is a general result from all three largest remnants of the inelastic 
collision (bouncing) simulations (17, 18, and 19). 



The same result is shown more quantitatively in the cumulative pressure plots in Fig. 10 which 
presents the results from inelastic collision simulations of 50-km targets. The interior materials 
(dotted line) have a steeper distribution than the exterior (dot-dashed), meaning that the interior 
experienced a more confined range of shock pressure compared to the exterior materials. The 
exterior material experienced both higher and lower shock pressures. However, although this result 
suggests that there may be a core of antipodal material that persists intact during the disruption 
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event, we find that this is not the case. The largest remnant is composed of material that originates 
from random radii in the target, hence any radial gradients in the original target are destroyed and 
not present in the largest post-collision remnant. 



3.4. Collision Models and Gravitational Reaccumulation 

Reaccumulation of remnants after a catastrophic disruption event involves slow to moderate 
velocity collisions and gravitational torques. The cloud of fragments is a mixture of dust, fractured 
fragments, competant fragments, and their gravitational aggregates. Each collision during reaccu- 
mulation is an impact calculation in its own right. This process must be computationally simplified. 



Benz and Asphaug (1999) and Melosh and Ryan (1997) used analytical methods to determine the 



mass of the largest post-collision remnant. In the simulations presented in this paper the largest 
remnant is determined directly from numerical integration in pkdgrav (Fig. [TJ. 

During reaccumulation, the physical processes that govern the outcome of collisions is assumed 
to be captured by the coefficient of restitution. We consider three different scenarios (as mentioned 



in {2.2), where the collision outcome is (i) perfectly inelastic (perfect merging, e n = 0), (ii) con- 
ditional merging (if Vi > V esc then inelastic bouncing with e n > else merging with e n = 0), and 
(iii) inelastic bouncing (e n > 0). In case one, merging replaces the two colliding particles with a 
single particle of the combined volume and average density, conserving angular momentum. Per- 
fect merging is the most computationally efficient, however, any particle that collides with another 
particle will merge regardless of the impact speed. Therefore, perfect merging may result in the 
largest post-collision remnant being unphysically massive. Conditional merging (case two) allows 
mergers only when the impact speed is less than the escape speed of the two particles, otherwise 
the particles rebound with energy loss determined by the coefficient of restitution. This is a more 
physical model; however, in both of the merging scenarios, information about the shape, surface, 
or internal structure of the remnants are not preserved. Inelastic collisions (case three) between 
particles is computationally expensive because the number of particles stays constant over the time 



of the simulation and the number of collisions that must be calculated is large (Leinhardt et al 



2000). Inelastic collisions preserves the shape and structural information about the remnants. 



Although all of these models are simplifications of the true reaccumulation process, perfect 
merging and inelastic bouncing represent the extremes, and they bound the true mass of the largest 
remnant. In Fig. [8] it seems that increasing the coefficient of restitution significantly decreases the 
critical catastrophic disruption energy. For example, under the same impact conditions, the largest 
remnant may be ~ 3 times larger using perfect merging versus inelastic bouncing (e.g., sim. 17a 
vs. 17b). The most striking difference between the two cases is the amount of most highly shocked 
material that reaccumulates into the largest remnant. 



However, closer examination of the end member cases reveals that the coefficient of restitution 
is not the controlling factor in the outcome of the pkdgrav calculation. In order to confirm the 
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robustness of the perfect merging case, all simulations were also run using the conditional merging 
criteria. The conditional merging results were nearly identical to the results of the perfect merging 
simulations, both in the mass of the largest remnants and the provenance of the material. We also 
tested the dependence of our bouncing simulation results on the value of the normal coefficient of 
restitution by running 3D simulations with two different values: e n = 0.5 (sim 17b, 18b) and 0.1 
(sim 17c, 18c). The different values for e n had no significant effect on the outcome of simulations 
using inelastic collisions. The lack of dependence on e n is due to the nature of the problem, where 
the majority of collisions occur between particles that are gravitationally bound to the largest 
remnant. 

Finally, we examined the criterion for reaching dynamical equilibrium in the reaccumulation 



of the largest remnant. As a comparison to the criteria described in § 2.2 we also used the more 



sophisticated algorithm in the companion code (Leinhardt and Richardson 2005). We calculated 



the amount of material bound to the largest remnant using the hierarchical system search feature of 
companion. The results are in excellent agreement, as we find that the mass of the largest remnant 
does not omit any significant mass from an extended bound system. 

We conclude that the results from the end member models (merging and bouncing) are com- 
putationally robust. Therefore, the differences in the merging versus bouncing results are due to 
the underlying physical differences in the models and not due to the value of a model specific 
parameter. 

A significant amount of the difference between the merging and bouncing results can be ex- 
plained as the result of runaway merging in both the perfect and conditional merging models. The 
first particles merge very early in the pkdgrav part of the simulation. The merged particles are 
replaced with a new particle that conserves volume at their center of mass. As a result of the mod- 
ified geometry upon merging, the new particle will very likely overlap neighboring particles, which 
in turn will most likely merge with the new particle. The process results in an artificial increase 
in the mass of the largest remnant. Note that all previous studies using the hybrid hydrocode-to- 
A^-body code technique have used the perfect or conditional merging criteria and suffer from the 
same artifact. 

On the other hand, the inelastic bouncing simulations probably underestimate the mass of the 
largest remnant. At the time of handoff between CTH and pkdgrav, particles at the base of the 
impact crater begin to bounce when in reality they would continue to move down into the target 
and shear (and not bounce at all). The abrupt transition from shearing flow to discrete, bouncing 
particles artificially disrupts a portion of the evolution of the problem. As a result, we believe 
that reality lies somewhere in between the merging and bouncing results. In future work, we will 
investigate these technical issues further by examining the modification in the material flow field 
from the handoff process. 
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Discussion 



Predicting the outcome of collisions in the outer solar system is challenging because there are 
several competing effects. Outer solar system bodies most likely have lower strength in comparison 
to asteroids of a similar size, yet they may be very porous which may help limit the damage from 



an impact (Asphaug et al. 1998 1 . Small outer solar system bodies should contain a high fraction 



of volatiles but this may be mixed with other materials, either intimately or segregated into layers 
(differentiated), which could have a complicated effect on the collision outcome. We have developed 
a new hybrid numerical technique that has the capability of tackling all of these complexities. In 
this work, we focus on nonporous bodies of varying shear strength and composition. 

Using our hybrid numerical technique, we derive new catastrophic disruption criteria for non- 



porous bodies of varying strength. Figure 11 presents a range of Q* D curves for outer solar system 
bodies, including both the strength and the gravity regimes. We focus on the mean impact velocity 
in the present-day Kuiper belt, about 1 km s _1 (Trujillo et al. 2001). The 3 km s _1 disruption 



curves for strong basalt and ice from Benz and Asphaug (1999) are presented as upper limits in 
strength and velocity for the Kuiper belt. As noted by Benz and Asphaug (1999), their 0.5 km s _1 



ice simulations produced results that were not consistent with the faster simulations and should be 
considered with caution until further work is done on the disruption of icy bodies. 

The critical disruption energy is a central component of models of the collisional evolution of 
a population of small bodies. Q* D curves are important because the transition from the strength to 
gravity dominated regimes defines the minimum required energy for disruption and the size of the 
weakest bodies. Thus, the amount of collisional grinding during the evolution of the solar system 



is sensitive to the disruption criteria (eg. Kenyon and Bromley 2004). In addition, the transition 
from strength to gravity produces observable kinks in the size distribution of small bodies, as found 



in the asteroid belt (e.g., O'Brien and Greenberg 2005) 



At present, the best Q* D curves for the outer solar system are still inadequate. Small bodies in 
the outer solar system are diverse, and we expect a wide range composition, bulk density, and shear 
strength. Here, we present our best recommendation for nonporous bodies. We will investigate 



porous and layered bodies in future work. Figure 11 shows a range of strength and composition 



from strong rock (Benz and Asphaug 1999) to weak ice targets. Here, we assume that changing 



the equation of state from solid rock to solid ice moves the Q* D curve down in the gravity regime 
by a factor of the bulk density, ~ 3 between the thick solid line to the thick dashed line. There 
are additional effects from changing equation of state, e.g., phase changes which modify the shock 
pressure decay profile; hence, the offset is not precisely the difference in density (as shown by 
Benz and Asphaug|[X9 99 ) . However, at the moment, these effects are within the numerical error in 
determining the mass of the largest remnant. We also focus on 90° impact events and note that 
changing the impact angle from 90° to and 45° is also within the range of errors in calculating Q* D . 
Note the small offset between the 90° points (open symbols) and ~ 45° thin lined curves in Fig. 11 
More oblique impacts (> 45°) have a more significant effect on Q* D ( Benz and Asphaug||1999 ). 



-22- 



In the strength regime, we recommend using laboratory experiments on ice targets (filled circles 



in Fig 


11 


from Giblin et al. 


2004| |Ryan et al. 1999| Kato et al. 


1992; 


Kato et al. 1 19951 |Arakawa 


1999 i 


\rakawa et al.||1995 


Cintala et al.||1985| Kawakami et al.||1983 


Lange and Ahrens||1987| Higa 



et al. 1996| |1998| ) and assume a strain rate dependent slope as described in Housen and Holsapple 



(1990|. Note that the strength regime of Benz and Asphaug (1999 1 's ice curves are much stronger 



than the laboratory data. The laboratory ice experiments span a large range in impact velocities (4 
to 3500 m s _1 ) and initial temperature (77 to 263 K). However, the results are not monotonic with 
impact velocity and the scatter likely represents differences in target preparation. There is limited 
data on weak rocks, but one study on the catastrophic disruption of pyrophyllite lies between the 



ice data and the strong rock data in the strength regime (Takagi et al. 1984 1 . 



In the gravity regime, the disruption criteria can be estimated for a range of impact conditions 
by considering the density of the material, p, the impact velocity, and coupling parameter, p. 



Housen and Holsapple (1990) derive 



Q 



D 



(1) 



where C is a constant that is fitted to our simulations and accounts for the effects of shear strength, p 
in kg m -3 , Rt in m, and Vi in m s . Based on our and previous numerical experiments, p ~ 0.45. 
Note that Eq. [T] applies only to hypervelocity collisions where a shock wave develops, which is 
expected for most collisions after the end of the accretionary phase in the solar system. We fit 
C = 2.6 ± 0.9 x 10~ 8 for weak targets (sims. 1, 8, 14) and C = 5.2 ± 1.7 x 10~ 8 for strong targets 
(sims. 2, 10, 16), using p = 2820 kg m -3 and the disruption thresholds and critical velocities in 
Table |3j Note that our weak rock results (thick solid line in Fig. 11) are also applicable to the 
asteroid belt. 

The intersection of the strength and gravity lines defines the transition region, which spans a 



few decades in size. Using the functional form of Benz and Asphaug (1999), 

Q* D = q b R a + pqgRP, 

our recommended parameters for ice targets (with a weak strength) at 1 km s _1 
m~ 3 (low temperature density), qi, = 20 J kg -1 , q g = 3.5 x 10 -6 J m 3 kg" 



(2) 

are: p = 930 kg 
0.4, and (3 = 1.3 

(thick dashed line in Fig. 11 ). eft is centered on the laboratory experiments, pq g is the value of the 



a 



gravity section of the curve at 1 m, a is taken to match Benz and Asphaug ( 1999 ) which is in good 



agreement with Housen and Holsapple (1999), and (3 is an average gravity regime exponent based 



on our and previous work. 

In this work, we also developed techniques to study the physical properties of the largest 
remnant, which will be used in future work to compare to observations of Kuiper belt objects. 
Insights into the composition and preservation of Kuiper belt objects may be acheived in future 
efforts that combine observations and modeling of collisional processes ( Leinhardt et al.|2008 ). Even 
though the range of shock pressures is low (e.g., Fig.[8]C), because of the expected physical properties 
of Kuiper belt objects, some shock modification may be observed. For example, thermodynamic 
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analyses of shock wave experiments on 100 K H2O ice indicate that peak shock pressures of 1.6 to 



4.1 GPa are required to produce incipient and complete melting, respectively (Stewart and Ahrens 



2005). The results from simulation group 19 indicate that catastrophic disruption events in the 



Kuiper belt reach conditions for melting ice. The effects of heating from the shock may produce 
changes in the composition of chemistry of outer solar system bodies. Because the reaccumulation 
process produces a final remnant with surface materials that have experienced a wide range of peak 
shock pressures, we suggest that surfaces of large rubble piles in the Kuiper belt may have complex, 
heterogeneous surface features (e.g., spectral or albedo variations). 

Our analyses indicate that outer solar system bodies may be disrupted more easily than indi- 
cated in previous simulations of catastrophic disruption. As a result, the collisional evolution of the 



outer solar system may have produced more mass loss from collisional grinding (e.g., Kenyon and 
Bromley||2004 1 . At the present time, the relative roles of dynamical and collisional sculpting of the 
outer solar system are poorly understood. The first attempts at combined dynamical-collisional 



models have just begun (Charnoz and Morbidelli 2007 and S. Kenyon, personal communication). 



Charnoz and Morbidelli (2007) found it difficult to reconcile observations of the size distributions 



of comets and Kuiper belt objects with their hybrid model of dynamical and collisional evolution 
of the outer solar system, and they suggested that outer solar system bodies could not have lower 



disruption energies than found by Benz and Asphaug (19991. In this work, we find that it is 



likely that outer solar system bodies require lower collision energies than assumed by Charnoz and 



Morbidelli (20071, illustrating the youthful stage of outer solar system collisional studies compared 
to studies of the asteroid belt. Recent work on the asteroid belt has illustrated how combined 
dynamical-collisional studies may be used to infer the complex evolution of the inner solar system 
dBottke et afl2005a|b| | . We are making progress toward similarly sophisticated studies of the outer 
solar system. 



5. Future Development 

The present work focuses on the gravity regime. In the future, we will further develop the 
CTH-pkdgrav technique to include modeling in the strength regime. The ROCK model is able to 
reproduce the fracture patterns in laboratory impact crater ing and disruption experiments; however, 
computational techniques are needed to define the size distribution of intact fragments. While the 
largest fragments may be directly resolved, unresolved fragments will need a statistical model to 
determine the size distribution (e.g., the Grady-Kipp model used in previous work). When the 
largest intact fragment is greater than a single pkdgrav particle, methods need to be developed 
to bind multiple pkdgrav particles together for the gravitational reaccumulation portion of the 
calculation. Modeling the strength regime will allow study of the transition from the strength to 
gravity regime and calculation of full Q* D curves for specific materials. 



6. Conclusions 



We present results from a series of catastrophic disruption collisions using a new hybrid 
hydrocode-to-iV-body numerical technique. This method has been developed to study the out- 
come of collisions between bodies with complex material properties, such as Kuiper belt objects. 
Our hybrid method allows us to directly calculate the mass of the largest gravitationally reaccu- 
mulated remnants. Hence, we also developed techniques to analyze the physical properties of the 
reaccumulated remnants from catastrophic collision events. Understanding the role of collisions in 
changing the composition and internal structure of Kuiper belt objects is important because KBOs 
are the best representatives of the planetesimals that accreted into the outer solar system planets. 
We find that material in the largest post-collision remnants experience a large range of peak shock 
pressures. We suggest that large reaccumulated bodies may have heterogeneous surface features 
and that interior and exterior properties may be different. 

In this work, we use rocky, asteroid-like targets to facilitate comparison between our results 
and that of previous work. We show that, when similar strength parameters are used, the derived 



catastrophic disruption criteria, Q* D , is in excellent agreement with previous work (Benz and As 



phaug 1999; Melosh and Ryan 1997). We also demonstrate that the shear strength of a body has 
a significant effect on the Q* D , even in the gravity regime, and we recommend new catastrophic 
disruption criteria for nonporous weak rocky and icy bodies. 
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7. Appendix 

Material parameters used with the ROCK model in CTH are presented in Table |AT| The ROCK 



model is based upon the strength model developed by Collins et al. (2004) for the SALEB hydrocode 
and implemented into CTH by Senft and Stewart (2007). Note that the total damage, D, is the sum 
of the shear and tensile damage and limited to a maximum value of 1. 
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Table 1: Principle parameters for different strength models. Simulations denoted by R use the 
ROCK model (otherwise GEO model). Note that Yq = 10 7 for strong ROCK models; full parame- 
ters for the ROCK model are given in the Appendix. 

Strength Y Y M Y T Sim. 

Model Pa Pa Pa Num. 

Hydro 9, 15 

Weak 10 6 10 7 10 5 1, 4R, 8, 14, 17, 18, 19 

Medium 10 7 10 s 10 6 13 

Strong 10 6 3.5 x 10 9 10 5 2, 6R, 7R, 10, 11R, 12R, 16 
StrongHT 10 6 3.5 x 10 9 10 7 3, 5R 
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Table 3: Catastrophic Disruption Threshold 



Sim. 


V cri t [km/s] 


[J/kg] 


1 


0.83 ± .1 


2.1(±0.5) x 10 2 


2 


1.4 ± .1 


6.3(±0.8) x 10 2 


3 


2.5 ± .3 


1.9(±0.4) x 10 3 


4 


0.97 ± .05 


2.9(±0.3) x 10 2 


5 


1.5 ± .1 


7.4(±0.9) x 10 2 


6 


1.9 ± .5 


1.1(±0.6) x 10 d 


8 


2.50 ±.04 


1.79(±0.06) x 10 3 


9t 


1.7 


8.5 x 10 2 


10 


4.8 ±.3 


6.6(±0.9) x 10 3 


11 


4.9 ±.3 


6.8(±0.8) x 10 3 


13 


3.0 ±.1 


2.5(±0.2) x 10 3 


14 


1.30 ±.07 


1.8(±0.2) x 10 4 


15t 


1.2 


1.5 x 10 4 


16 


1.9 ±.2 


4.0(±0.7) x 10 4 



+No error bars on simulation sets that contain only two simulations. All other simulation sets have three simulations. 

Table 4: Catastrophic Disruption Threshold 
Sim. R CTit [km] Q* D [J/kg] 
7 0.141 ±.003 1.6(±0.2) x 10 2 
12 1.13 ±.06 7(±1) x 10 3 



Table Al: ROCK model parameters. 



Parameter 


Weak 


Strong 


StrongHT 


Description 


Y (MPa) 


1 


10 


10 


shear strength at zero pressure 


Y M (GPa) 


0.01 


3.5 


3.5 


shear strength Von Mises limit 




1.2 


1.2 


1.2 


coefficient of internal friction (D=0) 


Y c (Pa) 











cohesion (D=l) 




0.6 


0.6 


0.6 


coefficient of friction (D=l) 


Y T (MPa) 


-0.1 


-0.1 


-10 


tensile strength 


T m (K) 


1392 


1392 


1392 


melting temperature at ambient pressure 




1.2 


1.2 


1.2 


thermal softening parameter 


Pbd (GPa) 


0.003 


2.94 


2.94 


brittle-ductile transition pressure 


Pb P (GPa) 


0.004 


4.11 


4.11 


brittle-plastic transition pressure 


V 


0.23 


0.23 


0.23 


Poisson's ratio 


Dso 











initial shear damage 


D T o 











initial tensile damage 
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Fig. 1. — Example of a 3D, 45 degree, impact between two basalt spheres (Rp = 14 km, Rt = 50 km, v. L 
= 1.8 km s _1 ) using the CTH-pkdgrav technique. Rows 1-2: CTH calculation (0-60 s) in cross-section along 
the y — plane. Rows 3-4: pkdgrav calculation (60 s - 140 h). Row 1 depicts the projectile (light grey), the 
target (dark grey) and the adaptive mesh. Row 2 shows the projectile and target (beige) and the pressure 
due to the impact (grey scale in Pa). 3A is a cross section (y = plane) of the target at 60 s after handoff 
from CTH to pkdgrav. Colors represent the peak pressure attained during the impact (logarithmic range of 
0.01 to 7 GPa). 3B-4C show the entire 3D object, zooming out as the disruption process proceeds. Scale bar 
is 50 km. The last frame shows only the largest reaccumulated, post-collision remnant which equilibrates to 
45% of the target mass in this simulation (sim 18b). 
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Fig. 2. — A. Initial tracer locations. B. Tracer locations and pkdgrav particles at handoff. Tracer and 
particles are color coded by peak pressure attained during the impact event. Note that projectile material 
lines the crater floor, shown by pkdgrav particles but not tracers. 
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Fig. 3. — Variation in mass of the largest post-collision remnant with handoff time from CTH to pkdgrav. 
R T = 10 km, R P = 0.83 km, and V l = 1.9 km s" 1 (sim. 8c). 
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Fig. 4. — A. Catastrophic disruption criteria (Q* D ) for basalt targets with varying strength models and 
impact velocities (Tables [3] & B. Q* D data adjusted to Vj = 3 km s _1 . Solid lines: fit to 3 km s _1 , 
90° impacts (slope = 1.26, |Benz and Asphaug||1999| |. Dashed lines: |Melosh and Ryan| ( [T997| Eq. 5 with 
Vi — 3 km s _1 (slope = 1.5). Data points, with ler errors, from this work: black dots (Sims. 1, 8, 14), 
triangle (4R) - weak targets (black dotted line, slope = 1.2); blue x (9, 15) - hydrodynamic targets; red 
dots (2, 10, 16) and open triangle (5R) - strong targets (red dotted line, slope = 1.3); green point (3), open 
triangles (6R, 11R), open squares (7R, 12R) - strongHT targets. 
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0.5 1 1.5 2 



Q/Q* D 

Fig. 5. — Mass of the largest post-collision remnant normalized to the mass of the target versus the specific 
impact energy normalized to the catastrophic disruption threshold. Each 2D simulation with constant mass 
ratio from Table [2] is represented by solid hexagon. Each simulation with varying mass ratio is represented 
by an open triangle. A least squares fit to the data produces a slope of —0.48 ± 0.02 (black line). The 
same linear relationship holds for strong to hydrodynamic basalt targets and impact velocities from 0.7 to 
5.6 km s . 
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Velocity at Handoff, (log 10 V(m/s)) Shock Pressure, (log 10 P(GPa)) 



Fig. 6. — A. Peak shock pressure as a function of distance (adjacent to the centerlinc) in the target for 
simulations with varying shear strengths but identical impact conditions. In all cases, an 0.83-km radius 
projectile hit a 10-km radius target at 3.7 km s _1 . B. Cumulative plot of the maximum shock pressure 
distribution (in mass fraction) in the total mass. C. Cumulative plot of the velocity distribution at handoff 
to pkdgrav. D. Cumulative maximum shock pressure distribution of material rcaccumulatcd into the largest 
remnant. 
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Shock Pressure (log 10 (GPa)) 

Fig. 7.— A. A 2D 90° impact into a 50 km target at V* = 1 km s" 1 , M lr /M T = 0.73 (sim. 14d). B. A 
3D 45° impact into the same target at Vi — 1.8 km s^ 1 , M; r /Mr = 0.74 (sim. 18a). Solid line is total mass, 
dashed line is the largest remnant, dotted line is second largest remnant. 
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Fig. 8. — Provenance plots for the largest remnant in 3D simulations of collisions into 50-km radius targets. 
Each sphere corresponds to a pkdgrav particle; colors represent peak shock pressure. Left column: perfect 
merging results. Right column: inelastic collision (bouncing) results. Simulation numbers and M; r /Mx are 
noted in parenthesis. 
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Fig. 9. — The largest remnant, Mi r /Mx = 0.21, from 3D simulation of 45° impact onto 50-km basalt 
target (sim. 17b). A. View of surface. B. View of interior, showing rear hemisphere cut along the Y = 
plane. Colors correspond to peak shock pressure. 
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Fig. 10. — Interior versus exterior shock pressure history of the largest remnant from 3D inelastic bouncing 
simulations of 50-km radii targets. Solid line: all material; dashed line: all of largest remnant; dash-dot: 
exterior of largest remnant (> 0.75i?; r ); dotted: interior of largest remnant (< 0.75i?/ r ). 
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Fig. 11. — Suggested range of Q* D curves for small outer solar system bodies, such as Kuiper belt objects 
Solid lines and squares 



Benz and Asphaug ( 1999 ) results: thin grey 
lines - 3 km s _1 angle averaged fits and open symbols - 90° simulations. Closed circles - 90° 



rock; dashed lines and circles - ice. 

laboratory 

catastrophic disruption data on ice (see text). Thick lines - recommended 1 km s _1 weak target catastrophic 
disruption criteria. 



